;+
; NAME:
;       ADD_GRS_INFO
;
; PURPOSE:
;       Add information about GRS-associated BGPS sources to the EMAF
;       structure.  Also, since this routine is called whenever a new
;       structure is generated, this routine adds flag for
;       star-formation activity from Dunham et al (2011, ApJ, 731,
;       90).
;
; CATEGORY:
;       EMAF
;
; CALLING SEQUENCE:
;       ADD_GRS_INFO [,ANGSEP=angsep]
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       ANGSEP -- [arcmin] Angular separation between BGPS peak flux
;                 and GRS cloud center to be used for assigning GRS
;                 association [Default = 10'].
;
; KEYWORD PARAMETERS:
;       SILENT -- Shut yer yap.
;
; OUTPUTS:
;       NONE  [Saved version of the emaf structure.]
;
; OPTIONAL OUTPUTS:
;       NONE
;
; MODIFICATION HISTORY:
;
;       Created:  10/29/12, TPEB -- Initial version, with code from
;                                   make_dist_comp_plot.pro.
;       Modified: 11/02/12, TPEB -- Forced sorting of Corrected GRS
;                                   Distance input to avoid buggies.
;       Modified: 11/13/12, TPEB -- Force GRS clouds within 0.5 kpc of
;                                   tangent distance to the Reid
;                                   tangent distance, too.
;       Modified: 11/29/12, TPEB -- Added code to insert
;                                   star-formation flag from Dunham et
;                                   al (2011, ApJ, 731, 90).
;       Modified: 12/03/12, TPEB -- Set default ANGSEP = 10'
;       Modified: 01/28/12, TPEB -- Reset GRS_TAN criteria; clouds
;                                   w/in 0.25 kpc of d_tan are set to
;                                   d_tan.  Tangent distances for C85
;                                   & R09 are nearly identical; this
;                                   should fix d_tan mismatches due to
;                                   vlsr issues.
;                                                     
;
;-

PRO ADD_GRS_INFO, ANGSEP=angsep, SILENT=silent
  
  ;; Parse
  IF n_elements(angsep) EQ 0 THEN angsep = 10.
  silent = KEYWORD_SET(silent)
  
  ;; Load in necessary files
  restore,'./ancillary/grs_kda_data.sav',/ver
  restore,'./emaf_paper/bgps_emaf.sav',/ver
  n = n_elements(emaf)
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF ~ exists THEN galactic_params 
  
  ;; Create arrays used for GRS identification
  dsep = fltarr(n)
  vsep = fltarr(n)
  grsi = lonarr(n)
  
  
  ;; Loop through EMAF to check for GRS co-incidence
  IF ~ silent THEN message,'Checking for GRS co-incidence...',/inf
  FOR i=0L, n-1 DO BEGIN
     gcirc, 2, emaf[i].l, emaf[i].b, grs.glon, grs.glat, dis
     mind = min(dis,gind)
     vsep[i] = abs(emaf[i].vlsr - grs[gind].vlsr)
     
     ;; A DSEP =-1 means there is no coincidence for this object
     dsep[i] = (vsep[i] LE grs[gind].delv*1.0) ? mind/60. : -1
     IF dsep[i] GT angsep THEN dsep[i] = -1
     grsi[i] = gind
     
  ENDFOR
  
  emaf.grs = dsep GE 0.
  IF ~ SILENT THEN message,'There are '+string(total(emaf.grs),format="(I0)")+$
                           ' GRS-associated sources.',/inf
  
  ;; Load emaf
  emaf.grs_clemdist = grs[grsi].dis
  emaf.grs_v        = grs[grsi].vlsr
  emaf.grs_l        = grs[grsi].glon
  emaf.grs_b        = grs[grsi].glat
  emaf.grs_21cm     = grs[grsi]._21cm_
  emaf.grs_angsep   = replicate(angsep,n)
  
  ;; Derive GRS distance based on Reid (2009) rotation curve
  
  ;; Read in GRS distance correction file
  readcol,'./emaf_paper/corrected_grs_distance.txt',bad_gd_cnum,new_gd,format='I,F',$
          COUNT=ncorrected
  IF ncorrected GT 0 THEN BEGIN
     ;; Sort GRS list
     gind        = SORT(bad_gd_cnum)
     bad_gd_cnum = bad_gd_cnum[gind]
     new_gd      = new_gd[gind]
     ;; Figure out where, etc...
     e_ind = WHERE_ARRAY(bad_gd_cnum,emaf.cnum,n_e)
     c_ind = WHERE_ARRAY(emaf.cnum,bad_gd_cnum,n_c)
     IF n_e NE n_c THEN message,'Drink some of that stout, chubby!'
     ;; IF n_e NE ncorrected THEN message,'Drink some of that porter, chubby!'
     IF n_e NE 0 THEN emaf[e_ind].grs_clemdist = new_gd[c_ind]
  ENDIF
  
  ;; Calculate tangent distances & how many at FAR KDIST
  grs_dtan = !MW.R0 * cos(emaf.grs_l * !dtor) / cos(emaf.grs_b * !dtor)/1.d3
  
  ;; If within 0.25 kpc of dtan, assign to dtan
  TAN_BUFF = 0.0
  grs_tan = abs(emaf.grs_clemdist-grs_dtan) LT TAN_BUFF AND emaf.grs
  grs_far = (emaf.grs_clemdist GE grs_dtan+TAN_BUFF) AND emaf.grs
  message,'N GRS at TANGENT d: '+string(long(total(grs_tan)),format="(I0)"),/inf
  message,'N GRS at far KDIST: '+string(long(total(grs_far)),format="(I0)"),/inf
  
  ;; Calculate R09 distance based on vlsr & grs_far
  FOR i=0L,n-1 DO BEGIN
     
     ;; Calculate grs_dist based on GRS KDA resolution
     IF emaf[i].grs THEN $
        IF ~ grs_tan[i] THEN $
           emaf[i].grs_dist = (grs_far[i]) ? $
                              KDIST(emaf[i].grs_l,emaf[i].grs_b,emaf[i].grs_v,$
                                    /FAR)/1.d3 : $
                              KDIST(emaf[i].grs_l,emaf[i].grs_b,emaf[i].grs_v,$
                                    /NEAR)/1.d3 $
        ELSE emaf[i].grs_dist = !MW.R0 * cos(emaf[i].l * !dtor) / $
                                cos(emaf[i].b * !dtor)/1.d3
  ENDFOR
  
  ;;==========================================================
  ;; Add star-formation flag from Dunham et al (2011, ApJ, 731, 90).
  
  ;; Find emaf elements in s
  s = read_bgps_csv()
  ind = where_array(emaf.cnum, s.cnum, ncom)
  
  d11 = read_mrt('./local/mrt/apj385558t2_mrt.txt')
  d11 = d11[ind]
  
  g3 = d11.ego GT 0 OR d11.rms GT 0
  g2 = ~g3 AND d11.r08 GT 0
  g1 = ~g3 AND ~g2 AND d11.glim GT 0
  g0 = ~g3 AND ~g2 AND ~g1
  
  emaf.sfa = 0*g0 + 1*g1 + 2*g2 + 3*g3
  print,m4_stat(emaf.sfa)
  save,emaf,filename='./emaf_paper/bgps_emaf.sav',/ver
  
END
